function out = f_mycenter(X,rad)
xi = 0;
for xx = -1:0.05:1
    for yy = -1:0.05:1
        xi=xi+1;
        count(xi,1) = 0;    
        count(xi,2) = xx;
        count(xi,3) = yy;
        for j=1:length(X)
            if sqrt( (xx-X(j,1)).^2 + (yy-X(j,2))^2)<rad
                count(xi,1)=count(xi,1)+1; % one point inside the circle                
            end
        end
    end
end
[Ysort,Ind]=sort(count(:,1),'descend');

% count

out.C1 = [count(Ind(1),2) count(Ind(1),3)];
i=2;flag = 0;
while i <=length(Ind) & flag == 0 % far from the 1st center
    i=i+1;
    % if sqrt( (count(Ind(1),2)-count(Ind(i),2)).^2 + (count(Ind(1),3)-count(Ind(i),3))^2)>rad*2
    if sqrt( (count(Ind(1),2)-count(Ind(i),2)).^2 + (count(Ind(1),3)-count(Ind(i),3))^2)>0.4
        flag = 1;
    end
end
out.C2 =  [count(Ind(i),2) count(Ind(i),3)];
  

% i=3;flag = 0;
% while i <=length(Ind) & flag == 0 % far from the 2nd center
%     i=i+1;
%     if (sqrt( (count(Ind(i),2)-out.C1(1)).^2 + (count(Ind(i),3)-out.C1(2))^2)>rad*2) & (sqrt( (count(Ind(i),2)-out.C2(1)).^2 + (count(Ind(i),3)-out.C2(2))^2)>rad*2)
%         flag = 1;
%     end
% end
% out.C3 =  [count(Ind(i),2) count(Ind(i),3)];


